function [ moments_y ] = recover_gauss_from_lognormal( moments_x )

    

% recovers the parameters of the Gaussian distribution used to generate a
% lognormal, based on the moments of the lognormal
% inputs are MEANS and VARIANCES
% outputs are also MEANS and VARIANCES
% in puts are mean_log_mu_i, var_log_mu_i, mean_log_v_i, var_log_v_i, log_rho_i
% for details see here: https://link.springer.com/chapter/10.1007/978-3-319-65112-5_20
% x is lognormal, y is normal

assert( numel(moments_x)==5, 'input must consist of 5 moments' ) 


%% 

[E_x1, V_x1, E_x2, V_x2, rho_x] = deal_cols( moments_x );

% means squared
Ex1_sq = E_x1^2;
Ex2_sq = E_x2^2;

% compute means of Gaussian
mu_y1 = log(  Ex1_sq / sqrt( Ex1_sq + V_x1 ) );
mu_y2 = log(  Ex2_sq / sqrt( Ex2_sq + V_x2 ) );

% compute variances of Gaussian
V_y1 = log(1+ V_x1/Ex1_sq );
V_y2 = log(1+ V_x2/Ex2_sq );

sd_x1 = sqrt(V_x1);
sd_x2 = sqrt(V_x2);
sd_y1 = sqrt(V_y1);
sd_y2 = sqrt(V_y2);

cov_x = rho_x * sd_x1 * sd_x2;

Cov_y =  log( 1 + cov_x/(E_x1*E_x2));
rho_y = Cov_y/(sd_y1*sd_y2);

moments_y = [ mu_y1, V_y1, mu_y2, V_y2, rho_y ];

end

